Two-dimensional Vesicle Dynamics under Shear Flow: Effect of Confinement 
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Dynamics of a single vesicle under shear flow between two parallel plates is studied in two- 
dimensions using lattice-Boltzmann simulations. We first present how we adapted the lattice- 
Boltzmann method to simulate vesicle dynamics, using an approach known from the immersed 
boundary method. The fluid flow is computed on an Eulerian regular fixed mesh while the location 
of the vesicle membrane is tracked by a Lagrangian moving mesh. As benchmarking tests, the 
known vesicle equilibrium shapes in a fluid at rest are found and the dynamical behavior of a vesicle 
under simple shear flow is being reproduced. Further, we focus on investigating the effect of the 
confinement on the dynamics, a question that has received little attention so far. In particular, we 
study how the vesicle steady inclination angle in the tank-treading regime depends on the degree 
of confinement. The influence of the confinement on the effective viscosity of the composite fluid is 
also analyzed. At a given reduced volume (the swelling degree) of a vesicle we find that both the 
inclination angle, and the membrane tank-treading velocity decrease with increasing confinement. 
At sufficiently large degree of confinement the tank-treading velocity exhibits a non-monotonous 
dependence on the reduced volume and the effective viscosity shows a nonlinear behavior. 

PACS numbers: 47.63.-b, 47.11.-j, 82.70.Uv 
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I. INTRODUCTION 

The study of blood flow at the microscale, i.e. the 
scale of blood corpuscules, is an important issue. In re- 
cent years this field has embraced several communities 
ranging from medical scientists to mathematicians. Clas- 
sical continuum approaches of blood flow, dating back to 
a century ago at least [l[, are based on several assump- 
tions and approximations that are both difficult to jus- 
tify or to validate. For example, in the microvasculature, 
where most of the blood flow resistance takes place, red 
blood cells (RBCs), which are by far the major compo- 
nent of blood, have a size which is of the same order 
as that of the blood vessel diameter. Thus, one expects 
that the discrete nature of blood should play a decisive 
role in microcirculation. A prominent example is the 
Fahraeus-Lindqvist effect: RBCs cross-streamline migra- 
tion towards the blood vessel center results in a dramatic 
collapse of blood viscosity, causing a reduction of blood 
flow resistance in the microvasculature. Even in larger 
blood vessels (e.g. veins, arteries) a satisfactory phe- 
nomenological continuum approach is lacking. One may 
thus hope that a constitutive law for blood will ultimately 
emerge from numerical simulations taking explicitly into 
account the blood elements. Still blood flow simulation 
is a challenging task since it requires solving for the dy- 
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namics of both the blood elements and the suspending 
fluid (plasma). 

Different numerical methods have been developed to 
study RBCs or their bimimetic counterparts (represented 
by vesicles and capsules), each having its own advantages 
and drawbacks. A widely used method is the boundary 
integral method which is based on the use of Green's 
function techniques 0. It has been successfully applied 
to vesicles The advantage of this method is the 

high precision. However, except for special geometries 
(e.g. unbounded fluid domain, semi- infinite domain), 
an appropriate Green's function is not available. This 
means that extra integrations over boundaries delimiting 
the fluid have to be performed, which increases the com- 
putational time significantly. In addition, this method 
is valid for Stokes flow only (no inertia) . Other classes 
of methods are phase-field 0-Q and level set approaches 
[l(J which can be applied both in the Stokes and Navier- 
Stokes regimes. Their advantage is the ability of han- 
dling, in principle, many particles by just specifying the 
initial condition (in any new run with different vesicles 
number only specifying initial conditions is required in 
principle!). However, these methods introduce a finite 
thickness of the membrane, which seems, up to now, to 
set a quite severe limitation regarding extraction of quan- 
titative data in the dynamical regimes. This requires a 
finite element technique with a grid refinement. Other 
types of methods consist of solving the fluid equations 
by adopting a "coarse-grained or mesoscopic" technique. 
Examples include the so-called multiparticle collision dy- 
namics (MPCD) or stochastic rotation dynamics (SRD) 
[Til [l2j . Its advantages are the relative ease of implc- 
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mentation and inherent thermal fluctuations which make 
the method very efficient if these are required. 

In this paper we apply an alternative mesoscopic 
method, namely the lattice-Boltzmann (LB) method. In 
the spirit of the LB method, a fluid is seen as a cluster 
of pseudo-fluid particles, that can collide with each other 
when they spread under the influence of external applied 
forces. Advantages of the LB method are its relative ease 
of implementation together with its versatile adaptabil- 
ity to quite arbitrary geometries. The LB method has 
been already adapted and used to perform simulations 
of deformable particles such as capsules [3l , vesicles and 
red blood cells [Til u nder flow. The main issue of the 
work presented in Iff is to accomplish simulations with 
a large number of particles while using a small number 
of nodes to reduce the computational time and the re- 
quired memory. This has been achieved by using ad hoc 
membrane forces that penalize any deviation from the 
equilibrium configuration. In the present paper we use 
the precise analytical expression of the local membrane 
force as it has been derived [l6j from the known Hel- 
frich bending energy (l7j . The perimeter conservation in 
our case is achieved by using a field of local Lagrangian 
multiplicators (equivalent to an effective tension). To 
accomplish the fluid-vesicle coupling we follow the same 
strategy used in [l4[ to simulate capsules dynamics. In 
flij the flow is computed by LB. The flow-structure two- 
way coupling is achieved using the Immersed Boundary 
Method (IBM). Although confining walls were considered 
in the above mentioned studies their effect on the dynam- 
ics was not studied. We believe that it is of interest to 
study the impact of the walls on the dynamics of vesi- 
cles, a question that to the best of our knowledge has not 
been treated in the literature so far for vesicles, capsules 
or red blood cells, but only for a droplet [5i| and a hard 
sphere [43j . 

Vesicles are closed lipid membranes encapsulating a 
fluid and are suspended in an aqueous solution. Their 
membrane is constituted of lipid molecules (also the ma- 
jor component of the RBC membrane) [18|. Each one 
has a hydrophilic head and a two hydrophobic tails. 
These molecules re-organize themselves if they are in con- 
tact with an aqueous solution, or properly speaking self- 
assemble, into a bilayer in which all the heads of the 
molecules are facing either the internal fluid or the ex- 
ternal one. Experimentally, vesicles with size of the order 
of 10/im - called Giant Unilamellar Vesicles (GUV) - can 
be easily prepared in the laboratory using, for example, 
the electro- formation technique [19|. Unlike RBCs, for 
vesicles we can vary their intrinsic characteristic param- 
eters (size, degree of deflation, nature of internal fluid, 
etc.). Despite the simplicity of their structure, vesicles 
have exhibited many features observed for red blood cells: 
equilibrium shapes l20l. ta nk-treading motion 0, |2ll . lat- 
eral migration 0, 122l.l23l| , or slipper-like shapes |24l . [25| . 
Capsules (a model system incorporating shear elasticity) 
have also revealed some common features with vesicles 
li. 



In the following sections we briefly introduce the for- 
mulation of the LB method for vesicles. We then study in 
two-dimensions (2D) the tank-treading motion of a single 
vesicle under shear flow between two parallel plates. Here 
we decided for 2D simulations since they are computa- 
tionally less demanding, but still capture all the relevant 
physics. We use large systems (in lattice units) because 
of the higher resolution required to extract the results 
shown below. Previous works done in 2D dealing with 
vesicles (also for red blood cells) have demonstrated that 
the dynamics in the third dimension is not relevant even 
in confined geometries [HI, [H[ . Vesicle dynamics un- 
der shear has been extensively studied in the literature. 
It is known that a vesicle placed in shear flow performs 
different kinds of motions depending on its degree of de- 
flation, the viscosity contrast between the internal and 
the external fluids and on the strength of the shear flow 
(see the phase-diagram in [11). When the viscosities of 
the internal and the external fluids are identical the vesi- 
cle performs a tank-treading motion. Its main long axis 
gets a steady inclination angle with respect to the flow 
direction while its membrane undergoes a tank-treading 
like motion. However, in the majority of the previous 
theoretical and numerical works the vesicle is placed in 
an infinite fluid (unbounded domain). This corresponds 
to the situation where the walls are too far from the vesi- 
cle to have any influence on its dynamics. For this reason 
here we study vesicle dynamics in a confined geometry. 
However, studying numerically the dynamics of vesicles 
in such conditions is a challenging problem from a com- 
putational point of view, especially in highly confined 
situations. We need to solve for the flow of the internal 
and the external fluids. The boundary separating the two 
fluids is also an unknown quantity since the membrane 
shape is not known a priori. 

Since the vesicle size (~ 10/im) is much larger than 
its membrane thickness (~ 5nm), mathematically we 
model the membrane as an interface with zero thickness. 
Tracking the motion of this freely moving interface under 
flow is not a simple task, especially when the membrane 
undergoes larger deformations due to hydrodynamical 
stresses. We need to label the interface by points which 
we track in time. Further, to take into account the de- 
formation an increased number of label points is required 
for the code to be stable and to capture deformation with 
good resolution. On the other hand, spatial derivatives 
on the membrane are needed to be evaluated at every 
time step to compute the membrane force. We need to 
evaluate the local curvature that is the fourth derivative 
of the vector position. Any formation of a highly buckled 
region in the membrane will introduce potential instabil- 
ity. Furthermore, the vesicle volume (the enclosed area 
in 2D) and its surface area (the perimeter in 2D) have to 
be kept conserved in time. At higher degrees of confine- 
ment possible undesirable contact between the membrane 
and the walls of the channel can be expected, and this 
is an additional difficulty to cope with. We do not use 
any ad hoc repulsive force from the wall, rather the non 



3 



contact is achieved via a proper handling of the viscous 
lubrication forces by the lattice Boltzmann method. 

We shall discuss how the vesicle-fluid coupling is ac- 
complished. For that purpose an approach known from 
the immersed boundary method [29[ is adopted. We 
present tests of the code by investigating vesicle equi- 
librium shapes in a fluid at rest. We then present simu- 
lation results regarding the steady inclination angle and 
the effective viscosity, as well as the tank-treading veloc- 
ity as functions of the reduced volume and the degree of 
confinement. 



II. THE LATTICE-BOLTZMANN METHOD 

The motion of the membrane can be induced by ex- 
erting an externally applied flow, and this is the physical 
situation we are interested in. In the present section we 
discuss how the fluid flow is solved for by using the LB 
method. In recent decades, the LB method has been 
introduced and widely used to simulate e.g. fluid flow 
in complex geometries (e.g. in porous media), multi- 
component and multi-phase flow (e.g droplets and bi- 
nary fluids) [3(| HH . Such popularity of the LB method 
among scientists and engineers has been gained thanks to 
its straightforward implementation and its local nature 
that allows for parallel programming. 

In the limit of small Mach Ma (ratio of the speed of a 
fluid particle in a medium to the speed of sound in that 
medium) and Knudsen Kn (ratio of the molecular mean 
free path to the macroscopic characteristic length scale) 
numbers the LB method is known to recover with good 
approximation the Navier-Stokes equations (30l.[3l|: 

p^^ + u-Vuj = -Vp + r/V 2 u + F, (1) 

Vu = 0, (2) 

governing the fluid flow of an imcompressible Newtonian 
fluid, p and rj are the mass density and the dynamic vis- 
cosity of the studied fluid, u and p are its velocity and 
pressure fields, and t is the time. F on the right-hand side 
is a bulk force (e.g. gravity) or the membrane forces as is 
the case for vesicles immersed in that fluid (see below). 
In the spirit of the LB method, a fluid is seen as a clus- 
ter of pseudo-fluid particles, that can collide with each 
other when they spread under the influence of external 
applied forces. In the LB context, not only the spatial 
position is discretized but also the velocity. This implies 
that every pseudo-fluid particle can move just along dis- 
crete directions with given discrete velocities. The main 
quantity associated with a pseudo-fluid particle is the dis- 
tribution function /i(r,i), with < fi < 1, which gives 
the probability of finding at time t the pseudo-fluid parti- 
cle at position r and having velocity c;, in the i-direction. 
There is no unique way in the choice of a lattice in the LB 
method. What matters is that the discretization has to 
fulfill the following constraints: mass conservation, mo- 
mentum conservation and isotropy of the fluid. Here, we 



adopt the so-called the D2Q9 lattice, where D2 is an ab- 
breviation for two-dimensional space while Q9 refers to 
the number of possible discrete velocity vectors [IH ■ 

The evolution in time of the distribution /; is governed 
by the LB equation: 

/i(r + Ci Ai, t + At) - /i(r, t) = At (A 5 + Fi) (i = 0...8), 

(3) 

where fi(r,t) is the old distribution of the pseudo-fluid 
particle when it was at position r at previous time t, 
/i(r + CiAt, t + At) is the new distribution of the same 
pseudo-fluid particle after it moved in the direction Cj to 
the new location r + CiAt during the elapsed time At, 
with At being the time step. The grid spacing is referred 
to by Ax. In this paper all units are given in lattice-units, 
where Ax = At = 1. The left-hand side of Eq. ([3]) alone 
represents the free propagation of the pseudo-fluid parti- 
cles without externally applied forces. In the right-hand 
side of Eq. F is any externally applied force and Aj 
is the collision operator. Here, we adopt the Bhatnagar- 
Gross-Krook (BGK) approximation which is given by 

A i = -^[/i(r,t)-/ i cq (r,t)]. (4) 

The BGK collision operator describes the relaxation of 
the distribution /i(r, t) towards an equilibrium distribu- 
tion fi q (r,t), with a relaxation time r. This relaxation 
time is set to 1 in this paper and related to the dynamic 
viscosity 77 via the relation: 

T] = vp = pc s -^ It- -I , (5) 

where c s = is the speed of sound for the D2Q9 

lattice. /j° q (r,t) is the equilibrium distribution obtained 
from an approximation of the Maxwell-Boltzmann distri- 
bution and is given by 




where uj-^ are weight factors. Wi equals 4/9 for the veloc- 
ity vector, 1/9 in the horizontal and vertical directions 
and 1/36 in the diagonal directions. The macroscopic 
quantities describing the flow are given by 
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p(r,t) = J2fM, (7) 

i=0 

for the local mass density, 

1 8 

u(r,f) = ^-7Y^/i(r,*)ci, (8) 
p{r ' l > i=0 

for the local fluid velocity and 

p(r,t) = p(r,t)c 2 s , (9) 
is the local fluid pressure. 
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The computational domain is a rectangular box, with 
length 2L and height 2W. We use x for the horizontal 
position of the box and y for the vertical position. Peri- 
odic boundary conditions are imposed on the right and 
on the left side of the box. To generate a shear flow, the 
upper and lower walls are displaced with the same veloc- 
ity u wa u but in opposite directions. To achieve this nu- 
merically, within the LB technique, the following bounce 
back boundary conditions are implemented on the two 
walls as [33j 

/_i(r,t + Af)=/ i (r,t) + 2^(u TOO «-ci). (10) 

Here, /_,: denotes the distribution function streaming in 
the opposite direction of i. In the absence of a vesicle, the 
flow relaxes towards a steady linear shear velocity profile 
of the form u°° = jyci, where 7 = u waU /W is the shear 
rate. 



III. FLUID- VESICLE INTERACTION 

We denote by f2 e xti ^int the fluid domains outside 
and inside the vesicle, respectively, and by T the vesi- 
cle boundary. The flow has to be computed considering 
boundary conditions on the membrane, which is a free 
moving interface. At the membrane T we require the 
continuity of the flow velocity 

u oxt (r m ) = u int (r m ) = v(r m ), with r m E T. (11) 

The ext and int suffixes are for the external and the in- 
ternal fluids, respectively, v is the velocity of any point 
r m belonging to the membrane. The continuity of the 
tangential velocities of the two fluids on each side of 
the membrane follows from the assumption of the no- 
slip boundary condition at the membrane. Continuity of 
the normal velocity is a consequence of mass conservation 
(integrating the incompressibility condition V • u = on 
a small volume straddling membrane and using the diver- 
gence theorem yields that condition). Continuity of the 
two fluid velocities with that of the membrane expresses 
the fact that the membrane is non permeable (normal 
component) and that we assume full adherence (tangen- 
tial component) (34[. Force balance (in the absence of 
inertia) implies that the net force acting on a membrane 
element is zero: 

(a cxt (r m ) - <r int (r m )) n = -f(r m ), with r m e T. (12) 

a is the hydrodynamical stress expressed by <7jj = 
rj(diUj + djUi) — pSij and n the unit vector normal to 
the membrane, pointing from the interior domain of the 
vesicle to the exterior one. f is the force exerted by the 
membrane on its surrounding fluid and its expression is 
given below. At sufficiently large distance from the vesi- 
cle membrane, the perturbation of the velocity field due 
to the membrane decays so that the fluid flow recovers 



its undisturbed pattern: 

u oxt (r) — > u°°(r), (13) 

|r— r m — >oo 

where r m <= T and r G Sl e xt- 

In what follows we show how these boundary condi- 
tions can be used to achieve the coupling between the 
fluid flow and the vesicle dynamics. In the present 
work the internal and the external fluid flows are com- 
puted by the LB technique. The velocity and the pres- 
sure fields are computed on an Eulerian regular fixed 
mesh, while the vesicle membrane is represented by a 
Lagrangian moving mesh immersed in the previous fluid 
mesh. The adopted method is the so-called immersed 
boundary method (IBM). This method was developed 
by Peskin to simulate blood flow in the heart |3]|. It is 
an adequate method to simulate deformable structures in 
flow (fluid-structure interaction) . For a review see for ex- 
ample Ref. [29| ■ Within the framework of this method an 
interface (separating two regions occupied by two distinct 
fluids) is discretized into points interconnected by clastic 
'springs' (as illustrated in Fig.[T]for the case of a vesicle). 
First, the fluid flow is computed in the whole computa- 
tional domain by ignoring the existence of the interface. 
Then, the interface is advected by the actual fluid veloc- 
ity, obtained from the Eulerian mesh by interpolation, as 
explained below. The fluid feels the existence of the vesi- 
cle due to singular point forces exerted by the interface 
nodes on their respective surrounding fluid nodes. This 
is achieved by linking the physical quantities computed 
on each mesh using a so-called discrete delta function 
suggested by Peskin (36[. The discrete delta function is 
defined as 

A «^( 1 + » s £)( 1 + c »2l:) < 14 > 

for ^ 2Ax and \y\ ^ 2Ax. In all other regions we set 
A(x) = 0, so that the function A has non zero values on 
a square. Here we choose 4Ax x 4Ax. The velocity at 
a given membrane node r m is evaluated by interpolating 
the velocities at its nearest fluid nodes r, n using the above 
A function so that we obtain 

v(r m ) =^A(r / -r m )u(r / ). (15) 
/ 

Here, u(r/) is obtained from the LB procedure. Deducing 
the velocity on the membrane nodes from the velocity of 
the fluid nodes is possible since we consider that the fluid 
velocity is continuous across the membrane and that the 
vesicle points are massless, behaving as tracer-like parti- 
cles, which do not disturb the flow at this stage. After 
evaluating every membrane node velocity we update its 
position using an Euler scheme 

r ro (t + Ai) =r m (i)+v(r m (i)), (16) 

and consequently the vesicle is advected and deformed. 
However, the vesicle membrane is not a passive interface. 
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FIG. 1: (Color online). Schematic view of a moving Lagrangian mesh representing a two-dimensional vesicle (where the 
membrane is represented by a contour) immersed in a fixed Eulerian mesh representing a fluid. 



It reacts back on the flow thanks to its restoring bending 
force 



membrane node. The force then takes the form 



f(r m ) = 



d 2 H 

ds 2 



2 



dC 

n+^t + KA {A-A )n, 

OS 



. (17) 

where H is the local membrane curvature, kb is the 
bending modulus, s is the arclength coordinate along the 
membrane (the contour in 2D), n and t are the normal 
and the tangential unit vectors, respectively. £ is a La- 
grange multiplier held that enforces local length conser- 
vation (the membrane is a one dimensional incompress- 
ible fluid). A detailed derivation of this force can be 
found in [16| ■ The additional last term in Eq. (fl7|) is in- 
troduced in order to enforce area conservation, because 
numerically a slight variation is observed (see Ref. 
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Aq is the initial reference enclosed area of the vesicle, A 
is its actual area and ka a parameter that is chosen in 
a such a way to keep the vesicle area conserved. This 
conservation constraint can be improved also by increas- 
ing the resolution. The membrane force has non zero 
value only on the membrane and should vanish elsewhere. 
More precisely, a given fluid point rj is subject to the 
force 



F(r/) 



f(r m )<5(r/ - r m )ds(r m ), with r m e T, 



(18) 

where ds is the distance between two adjacent mem- 
brane points. However, since the membrane is discretized 
and thus presented by a cluster of points, this integral 
is rather a sum of the singular forces localized on the 
membrane nodes. In addition, writing the force felt by a 
fluid node in terms of an ordinary Dirac delta function 
is not adapted here, since the membrane nodes can be 
off-lattice and do not necessarily coincide with the fluid 
lattice nodes. The Dirac delta function in Eq. (fT5)) is 
replaced by the A function suggested above which has 
a peak on the membrane node and decays at a distance 
equal to twice the lattice spacing after which it vanishes 
|36|. In this way the membrane force has a non zero 
value in a squared area of 4Ax x 4Acc centered on the 



F(r/) 



n 

E 

m— 1 



f (r. m )A(r/ - r m ). 



(19) 



Here, n is the number of membrane nodes. The vesicle 
membrane finds itself, on the one hand, advected by its 
surrounding fluid and on the other hand it exerts a force 
in response to the applied hydrodynamical stresses, caus- 
ing thereby a disturbance and modification of the fluid 
flow. 



IV. SIMULATION RESULTS AND DISCUSSION 

A. Dimensionless numbers 

The fluid flow and vesicle dynamics are controlled by 
the following dimensionless parameters: 

• The Reynolds number 

fnRl 



Re 



(20) 



is associated to the applied shear flow and mea- 
sures the importance of the inertial forces versus 
the viscous ones. Rq is the effective vesicle radius. 
In 2D, R can be deduced from the vesicle perime- 
ter i?o = -P/2-7T. In our simulations we use small 
enough values for Re (see below). 



The capillary number 



Ca = 



(21) 



represents the ratio between the shear time (I/7) 
and the characteristic time (wRq/kb) needed by a 
vesicle to relax towards its equilibrium shape af- 
ter flow cessation. This parameter controls the de- 
formability of the vesicle under flow. Larger Ca 
lead to a larger deformability. Below we use the 
value Ca — 1 which corresponds to the intermedi- 
ate regime. 
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FIG. 2: (Color online). Computed equilibrium shapes for vesicles having the same perimeter, but different reduced volumes 
a. Red colored solid lines are shapes computed by the LB method. For comparison purpose and for validation, we plot also 
equilibrium shapes computed by the boundary integral method (4Tj] (black dashed line). 



• The reduced volume (the swelling degree) quantifies 
how much a vesicle is swollen. In two dimensions 
it is given by 



A 
Ac~ 



p2 ' 



(22) 



where Ac is the area of a circle having the same 
perimeter P as the vesicle, a is unity for a circular 
vesicle (a maximally swollen vesicle) and less than 
unity for a deflated one < a < 1. 

• The viscosity ratio between the internal and exter- 
nal fluids is given by 



A 



Vint 
"text 



(23) 



In this paper, however, this ratio is taken to be 
unity. For this value, a vesicle is expected to un- 
dergo tank-treading motion 0, H3, HI] . 



The tension 



iber 



Ca s 



k p 



(24) 



is the ratio between the spring relaxation time (re- 
call that Kp is the spring constant) and the shear 
time (1/7). This number controls the inexten- 
sibility of the membrane under flow. To ensure 
the vesicle perimeter conservation constraint we set 
Ca s significantly small as compared to Ca (be- 
low we set Ca s = 1.05 x 10~ 5 ). For the simu- 
lations, we tune Kp until we get very negligible 
variations of the perimeter P. Kp is related to £ 
(the Lagrange multiplier) via the formula £(s, t) = 
Kp [ds(s,t) — ds(s, 0)1, where ds(s, 0) is the initial 
reference value [la, |34j . 

The degree of confinement is given by the ratio 
of the vesicle's effective radius to the channel half 
height, 



B. Computed equilibrium shapes 

Finding the vesicle equilibrium shapes constitutes one 
of the benchmarking tests we use to validate our code. 
In contrast to a droplet, which adopts spontaneously a 
spherical equilibrium shape, vesicles can adopt different 
kinds of non-spherical shapes. In two-dimensions, a vesi- 
cle gets a circular equilibrium only for a = 1. Usually 
the equilibrium shapes are obtained by minimizing the 
Helfrich bending energy [ItJ 



E = ^-jj2Hfds, 



(26) 



x = Ro/w. 



(25) 



subject to the two constraints of vesicle area A and 
perimeter P conservation (in 2D). The only parameter 
controlling the shape of a vesicle, in the absence of an 
external applied flow and in unbounded domain, is its re- 
duced volume a [20j . An alternative to energy minimiza- 
tion is to set the flow to zero (Re = 0) and let the vesicle 
relax to its terminal shape. Technically, we place initially 
a vesicle with some shape (here an ellipse) in a fluid at 
rest (no shear flow). Then, the membrane starts to de- 
form in order to relax towards the shape that minimizes 
its energy (Eq. l26|) . During this transition the membrane 
induces some weak fluid flow, inside and outside the vesi- 
cle. This flow stops once the vesicle gets its equilibrium 
shape. Figure [5] shows the computed equilibrium shapes 
for five vesicles having different values of the reduced 
volume a = 0.6, 0.7, 0.8, 0.9, and 1. The five vesicles 
have the same perimeter. Varying the reduced volume is 
achieved only by varying the vesicle area. It is somehow 
like swelling or deflating these vesicles to get different 
equilibrium shapes. To perform simulations we used the 
physical parameters Re = Ca = (fluid at rest), Rq = 20 
(to achieve a sufficient resolution at the scale of the LB 
grid) and x = Ro/W = 0.1. We have set n = 100, a value 
for which the code is stable. Significantly larger values 
of n may cause instability. From this point of view, the 
LB method differs from the other conventional numerical 
schemes, for instance the boundary integral or the finite 
difference/clement methods. In those methods, higher 
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resolution and higher stability is achieved by increasing 
(without limit) the number of discretization points. In 
contrast, with the LB method an increase of n induces 
higher resolution, but care should be taken not to exceed 
some given threshold value, beyond which the code desta- 
bilizes [ill] . Therefore, in all our simulations we have kept 
a sufficiently small enough number of membrane nodes 
per lattice grid (by keeping the distance between two ad- 
jacent membrane nodes ds close to 1). Within the LB 
method the velocity has to be kept small enough (in our 
case we choose the limit of 0.1) in order to have a suf- 
ficiently low Mach number and to ascertain the limit of 
neglectable fluid compressibility. The other parameters 
are chosen as follows. We have set L = 200, so that the 
flow perturbation due to the presence of the vesicle is 
negligible at the computational domain boundaries where 
periodic boundary conditions arc imposed, ka — 0.01 to 
fulfil a precise enough conservation of the vesicle enclosed 
area (we measure a variation of the order of 0.00015%) 
and we set up = 12 to keep the perimeter conserved as 
well (variation of 0.00125% is measured). The obtained 
shape for every given reduced volume is compared with 
its corresponding shape obtained by the boundary inte- 
gral method, the black dashed lines in Fig. [2] (the same 
method as used in Refs. 

[M S3)- For a given reduced 
volume, the computed equilibrium shapes obtained by 
both numerical methods are indistinguishable, especially 
at higher values of the reduced volume. In Fig. [2] we can 
see that for a reduced volume of 0.6, the vesicle assumes 
a biconcave shape, as it is typical for healthy red blood 
cells. 



(a) (b) 




-150 -75 75 150 -40 -20 20 40 

(c) x (d) 




Reduced volume Reduced volume 



FIG. 3: (Color online). Physical quantities associated to the 
tank-treading motion of vesicles (with Rq — 30) under shear 
flow (with Re = 9.45 x 1CT 2 and Ca = 1) in a confined 
channel: (a) streamlines pattern inside and outside a vesicle 
(a = 0.9) performing a tank-treading motion in a confined 
channel, (b) steady shapes for different values of the reduced 
volumes, (c) inclination angle versus the vesicle reduced vol- 
ume for two degrees of confinement \ ~ Ro/W = 0.40 and 
0.81, (d) membrane tank-treading velocity (scaled by jRo/2, 
the rotational velocity of a circular vesicle under shear flow 
in unbounded geometry) versus the reduced volume. 



C. Tank-treading under shear flow 

In the present section, we treat the effect of confine- 
ment on the dynamics of a tank-treading vesicle. First we 
study how the physical quantities, associated to the tank- 
treading regime, vary with the reduced volume. Then, for 
a given reduced volume, we analyze the effect of confine- 
ment on dynamics and rheology. We consider a single 
vesicle placed in a fluid subject to a simple shear. Here, 
we set Rq = 30 in order to achieve a high enough res- 
olution. For i? = 30, our explorations led us to the 
conclusion that n = 150 is a good compromise between 
numerical stability and resolution. For this value of n the 
code is stable even at higher degree of confinement. This 
also allows us to keep a sufficient number of fluid nodes 
between the wall and the membrane, a precision required 
in more confined situations. The length of the simulation 
box is set to L = 600, chosen as to minimize perturba- 
tions by the vesicle at the edge of the simulation box, 
where periodic boundary conditions are imposed. Under 
such conditions and in the absence of a viscosity con- 
trast (A = 1) a vesicle performs a tank-treading motion 
@, I371.I31I . It deforms until reaching a steady fixed shape 
with its main axis assuming a steady inclination angle 
with respect to the flow direction. The membrane un- 



dergoes a tank-treading like motion and so it generates 
a rotational flow of the internal enclosed fluid. 



D. Effect of the reduced volume 

Figure [3] shows different physical quantities measured 
in the tank-treading regime. In Fig. [3^ we show a vesi- 
cle performing tank-treading motion in a confined chan- 
nel. The vesicle assumes a steady inclination angle (the 
red solid line). The streamlines inside and outside the 
vesicle (the gray solid lines) show that the internal fluid 
undergoes a rotational flow, induced by the membrane 
tank-treading. The external fluid exhibits recirculations 
at the rear (the left side of the figure) and at the front 
(the right side of the figure) of the vesicle. Such recir- 
culations do not take place in the unbounded geometry 
Q. For a tank-treading vesicle in unbounded geometry 
(or at a sufficiently weak confinement), the external fluid 
lines are curved around the vesicle without being sepa- 
rated. In Fig. [3Jt the external fluid lines are separated 
into two portions before approaching the vesicle at two 
saddle points (located close to the channel centerline at 
the back and at the front of the vesicle): one portion con- 
tinues its flow (through the region between the wall and 
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FIG. 4: (Color online). Induced pressure field (with the grey scale) and flow streamlines (the gray solid lines in the right figure), 
inside and outside the vesicle, for different values of the reduced volume a = 0.6, 0.8 and 1. Re = 9.45 x 10~ 2 , Ca = 1 and 
X = 0.81. The black solid lines in the left figures and the red solid lines in the right figures represent the vesicle membrane. 
In the right figures, the regions with black color correspond to higher pressure while the white regions correspond to lower 
pressure. 



the membrane) and passes the vesicle, while the other 
portion is reflected back by the vesicle. Such flow re- 
circulations are also observed for confined rotating rigid 
spheres [4l[ and rigid ellipsoids [42| . For the same de- 
gree of confinement (\ = 0.4), in Fig. we varied the 
reduced volume of the vesicle. In Fig. |3Jd we report the 
steady state shapes obtained for different values of the 
vesicle reduced volume. All the vesicles in Fig. [3Jd have 
been initialized with a zero inclination angle. Figure [3b 
shows the steady inclination angle as a function of the re- 
duced volume (for two confinements: x = 0.4 and 0.81). 
The steady inclination angle increases monotonically (for 
both values of x) with increasing the reduced volume, un- 
til approaching 45 degrees in the limiting case of circular 
vesicles. The same qualitative tendency is observed in 
the unbounded geometry [H, H3, S3| ■ 

Figure [3Ji shows the behavior of the tank-treading ve- 
locity normalized by jRq/2, which is the tank-treading 
velocity of a circular vesicle [Hj|. For \ = 0.4, the tank- 
treading velocity increases monotonically with increasing 
the reduced volume, as observed for the unbounded ge- 
ometry [3l [37l |40|. However, for higher confinement, for 
example \ — 0.81, the tank-treading velocity does not 
vary anymore in a monotonous way. It has a maximum 
around a = 0.85 before it decreases at larger a. This be- 
havior can be explained by the fact that at higher degree 
of confinement, the amount of the external fluid which 
is able to flow from one side (the left) to the other side 
(the right) of the channel by crossing the narrow region 
between the wall and the membrane becomes smaller and 



smaller when increasing the reduced volume. At higher 
reduced volumes the inclination angle increases and so 
the membrane comes in closer proximity of the wall, see 
Fig. |4j Therefore, the external fluid flow does not par- 
ticipate fully to generate the tank-treading motion of 
the vesicle. This is also corroborated by the fact that 
the external fluid undergoes recirculation (see Fig. [5Jd 
and |4]) meaning that part of the fluid is reflected back- 
wards when approaching the vesicle. For a circular vesi- 
cle (a = 1) the amount of the reflected fluid is larger 
compared to the amount crossing the narrower region 
between the membrane and the walls. The induced 
pressure field shown in the three left panels in Fig. @] is 
significantly affected by increasing a. For a = 1, a sig- 
nificant pressure gradient is observed at the inlet and the 
outlet of the narrower gap between the membrane and the 
wall. Such pressure drop along this narrower region gen- 
erates an almost parabolic velocity profile, for the case of 
a = 1, as is shown in Fig. [5] In the same figure, for com- 
parison purpose, we report the disturbed velocity profile 
for other different values of the reduced volume (these 
profiles are taken at x = 0). The disturbance is maximal 
for a = 1. 

The bounce back boundary condition of Ladd [H[ al- 
lows to measure directly the hydrodynamical stress field 
a xy exerted by the fluid upon the channel walls. Fig- 
ures^ and[Bj3 show the measured hydrodynamical stress 
for two degrees of confinement, \ = 0.4 and 0.81. The 
effective viscosity 77 c ff can be extracted from the hydro- 



9 



dynamical stress using the formula 



where (u xy ) refers to the volume average of the stress ten- 
sor. The stress a xy has been averaged along the bottom 
wall. As shown in Figs. [BJi and [Ht> the hydrodynami- 
cal stress on the bottom wall exhibits important varia- 
tions for a larger reduced volume a. For a = 1, the 
stress is symmetrical with respect to the vertical axis at 
x = (which is perpendicular to the walls and passing 
through the center of mass of the vesicle). This symme- 
try is also observed for confined rigid spheres [43| , and is 
a consequence of the symmetry of the Stokes equations 
upon time reversal. Actually, our simulation contains 
a small amount of inertia, but so small that the asym- 
metry is difficult to identify on the figure. For deflated 
vesicles (a / 1) the stress curve has two unequal max- 
ima and one minimum. The flow deforms the vesicle and 
breaks the up-stream/down-stream symmetry (see Fig. [4] 
for a = 0.6 and a = 0.8). In other words, the Stokes 
reversibility is broken by the shape deformation. The 
values of these maxima and minimum significantly devi- 
ate from the corresponding value 777 in the absence of 
the vesicle (presented by the horizontal grey solid line in 
Fig. IB]) upon increasing the reduced volume. By compar- 
ing Figs.|B^, and[S]3, one notices that the stress is impor- 
tant for higher degrees of confinement, as expected. Sur- 
prisingly, at higher Rq/W we observe regions with nega- 
tive hydrodynamical stress. We believe that this results 
from a subtle effect due to fluid recirculation around the 
vesicle. However, a clear explanation of this phenomenon 




Flow Velocity 



FIG. 5: (Color online). Disturbed flow velocity profile mea- 
sured at x — for different values of the vesicle reduced vol- 
ume. The black dashed line is the undisturbed applied shear 
flow profile u x — 73/ in the absence of the vesicle. The pink 
solid line corresponds to the location of the membrane for 
a = 1. 



is at present not available. 

Figure [BJ: shows the behavior of the effective viscosity 
for different values of the vesicle reduced volume. The 
viscosity increases when increasing the reduced volume. 
The same tendency was observed for a vesicle placed in 
unbounded domain [Hj]. This result is explained as fol- 
lows: for a given confinement, the increase of the reduced 
volume implies a larger cross section of the vesicle in the 
channel (because of the large increase in the inclination 
angle), and this opposes more resistance to the fluid flow. 



E. Effect of confinement 

Here we set (a = 0.9) and vary only the width (2W) 
of the channel in order to study the effect of confine- 
ment. All the other physical and numerical parameters 
are kept identical to those of the previous section. All 
simulations in Fig. [Jj arc performed with \ varying from 
0.4 to 0.81. For this set of parameters, the code is still 
stable and guarantees a quite satisfactory conservation 
of the vesicle area and perimeter (AA/A ~ 0.00015% 
and AP/Pq ~ 0.013%). Note that in Figure [7J we have 
kept the same shear rate 7 = 1.75 x 10 -5 . Figure [7^, 
shows a decrease of the inclination angle upon increasing 
confinement. Under confinement the angle saturates at 
smaller values as compared to the corresponding one in 
the unbounded flow. The wall acts via a hydrodynami- 
cal repulsive force (a viscous force) tending to push, so to 
speak, the orientation angle back to zero in order to align 
the vesicle with the flow. Such a decrease of the inclina- 
tion angle was also reported by Janssen and Anderson 
[45[ | for a confined droplet under shear flow. 

By decreasing confinement further, we expect to reach 
the unbounded geometry limit (x = 0). We did not at- 
tempt to study this asymptotic limit. For x < 0.4 and 
within the present resolution (i?o = 30) a significant in- 
crease of L is required in order to avoid unphysical effects 
induced by the periodic boundary conditions. This task 
requires a very high amount of computational time. 

Figure [JJa shows a decrease of the membrane tank- 
treading velocity (scaled by jRo/2) versus confinement. 
At high enough confinement, the external fluid does not 
anymore participate wholly in membrane tank-treading. 
The fluid is partially reflected backwards when approach- 
ing the vesicle, an effect that increases with confinement 
(see Fig. [8]) . For rigid spheres a similar decrease of the 
rotation velocity is also observed when increasing con- 
finement [43l ]. Varying confinement affects also the way 
the membrane and the wall interact. Figurc[7J; shows the 
hydrodynamical stress exerted by the fluid on the bottom 
wall. The horizonal black solid line is the stress measured 
in the absence of the vesicle (iff). At distances far from 
the location of the vesicle (on the extreme right and left 
sides of the figure) the stress measured for all degrees of 
confinement matches with the stress in the absence of the 
vesicle. In the vicinity of the vesicle, around x = 0, we 
see deviation of the stress from the value 777. Such devi- 
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FIG. 6: (Color online). The hydrodynamical stress exerted by the suspending fluid upon the bottom wall for two degrees of 
confinement, \ = 0.4 (a) and 0.81 (b). The grey solid line is the stress calcuted analytically in the absence of the vesicle 777. 
(c) The deduced effective viscosity of the fluid, in the presence of the vesicle, for both degrees of confinement. 
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FIG. 7: (Color online). Variation of the measured physical 
quantities associated to a vesicle performing tank-treading 
motion in confined geometries when varying the degree of 
confinement, (a) steady inclination angle, (b) the mem- 
brane tank-treading velocity (scaled by "/Ro/2), (c) the hy- 
drodynamical stress field applied on the bottom wall and (d) 
the effective viscosity. Parameters are a — 0.9, Ro = 30, 
Re = 9.45 x 10~ 2 and Ca = 1. 



ations become larger and larger when increasing confine- 
ment. Again, as in the previous section, we observe stress 
with negative values at higher confinement. The effective 
viscosity is extracted from the stress (using Eq. [27j) and 
the results are shown in Fig. [TJi. The effective viscosity 
is found to significantly increase non-linearly with con- 
finement. 



In order to gain further insight we represent the pres- 
sure field and the streamlines (Fig. [5]), for three different 
confinements, x = 0.81, 0.6 and 0.4. Confining further 
the vesicle results in an increase of the pressure inside 
the vesicle, entailing a higher pressure gradient along the 
fluid layer located between the membrane and the wall. 
The amount of the fluid crossing this region decreases 
also when increasing confinement. The streamlines pat- 
tern shows that the recirculation becomes important at 
higher degree of confinement. Their two focal points 
move closer and closer to the vesicle when increasing con- 
finement. A closer inspection of the pressure field and the 
streamlines (for a given degree of confinement \) reveals 
interesting dynamics occurring in the narrow region be- 
tween the vesicle and the wall. When the external fluid 
approaches the vesicle it splits into two parts. One part 
is reflected by the vesicle and pushed backward without, 
passing the vesicle. The other part continues its flow and 
crosses the narrower gap formed between the wall and 
the membrane. At the inlet, of this gap, the pressure is 
significantly higher, resulting in a slowing down of the 
fluid (the streamlines are separated). Once the fluid en- 
ters this region its velocity is amplified (the streamlines 
come closer), under the action of the pressure gradient 
along the gap, until it exits that region. At the outlet, 
the pressure drops to a lower value and the fluid is slowed 
down again (the streamlines separate again). 

Finally, some general comments are worth mentioning. 
The fluid motion in the narrower gap, formed between 
the vesicle membrane and the wall, is induced under the 
action of three mechanisms: 1 - the membrane force, 2 
- the shear flow and 3 - the pressure-gradient along this 
region. The third mechanism dominates at high confine- 
ment. In that regime the fluid (in the gap) is subject to 
the sum of forces induced by the above three mechanisms. 
Within the present method, this sum must not exceed 
some given threshold otherwise the code becomes unsta- 
ble (due to higher flow velocities). There is also another 
technical detail that becomes problematic in situations 
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FIG. 8: (Color online). Induced pressure field and flow streamlines, indside and outside the vesicle, for different values of the 
degree of confinement \ = 0.81, 0.6 and 0.4. Parameters are a = 0.9, R = 30, Re = 9.45 x 10~ 2 and Ca = 1. 



of higher degrees of confinement. We assumed that the 
membrane has a zero thickness. However, by using the 
expression Eq. [TJ] the membrane force is distributed on 
fluid nodes located at distances of roughly 4Ax from the 
membrane. The membrane acquires a non-zero effective 
thickness. In more confined situations we need to leave 
at least 4 fluid nodes in the gap between the membrane 
and the wall, otherwise the dynamics of the vesicle suffer 
from numerical artifacts. For example, a leakage of the 
internal fluid is observed, and the tank-treading velocity 
exhibited a non uniform behavior along the membrane. 
To overcome all these problems we had to increase the 
resolution. For the resolution of Rq = 30 used in this 
section, the upper limit of the confinement we were able 
to reproduce without any apparent problem is \ = 0.81. 



V. CONCLUSION 

We have studied the effect of confinement between two 
parallel walls on vesicle dynamics under shear flow. We 
limited ourselves to the case of having the same fluid 
inside and outside the vesicle. In such a situation the 
vesicle performs tank-treading motion. We developed 
a latticc-Boltzmann method to perform two-dimensional 
simulations. The coupling between fluid flow and vesi- 
cle dynamics was adopted from the immersed boundary 



method. Unlike previous works, we have introduced the 
membrane force by using its analytical expression as a 
function of the mean curvature and its derivative. The 
vesicle enclosed area and its perimeter are kept conserved 
in our method. We first computed the known vesicle 
equilibrium shapes for different values of the swelling de- 
gree in order to validate our code. The obtained shapes 
match perfectly the ones computed by the boundary in- 
tegral method. As a second step, wc studied the case 
of a vesicle placed in a domain bounded by two paral- 
lel walls. We induced the shear flow by moving these 
two walls in opposite directions. We found that both the 
vesicle inclination angle, with respect to the flow, and 
its membrane tank-treading velocity decrease when in- 
creasing the degree of confinement. Moreover, since at 
sufficiently large degree of confinement the vesicle mem- 
brane comes close to the wall so that just a very narrow 
region is left for the external fluid to flow. Therefore, the 
vesicle acts as an obstacle and thus the effective viscosity 
increases dramatically when increasing confinement. At 
a given degree of confinement, we varied the swelling de- 
gree. We observed the same qualitative tendency as for 
the unbounded geometry for the behavior of the angle, 
tank-treading velocity and viscosity as a function of the 
swelling degree. However at higher degree of confinement 
even the angle still shows an increase with increasing the 
swelling degree, the measured values are lower. The tank- 
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treading velocity does not increase monotonically with 
the swelling degree. It exhibits a maximum value before 
getting to lower values in the limit of circular vesicles. 
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